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ABSTRACT 

We present a novel method to determine eccentricity constraints of extrasolar 
planets in systems with multiple transiting planets through photometry alone. Our 
method is highly model independent, making no assumptions about the stellar pa- 
rameters and requiring no radial velocity, transit timing or occultation events. Our 
technique exploits the fact the light curve derived stellar density must be the same 
for all planets transiting a common star. Assuming a circular orbit, the derived stellar 
density departs from the true value by a predictable factor, 5*, which contains in- 
formation on the eccentricity of the system. By comparing multiple stellar densities, 
any differences must be due to eccentricity and thus meaningful constraints can be 
placed in the absence of any other information. The technique, dubbed "Multibody 
Asterodensity Profiling" (MAP), is a new observable which can be used alone or in 
combination with other observables, such as radial velocities and transit timing vari- 
ations. An eccentricity prior can also be included as desired. MAP is most sensitive 
to the minimum pair-combined eccentricity e.g. (ei -I- e2)min- Individual eccentricity 
constraints are less stringent but an empirical eccentricity posterior is always derivable 
and is freely available from transit photometry alone. 

We present a description of our technique using both analytic and numerical im- 
plementations, followed by two example analyses on synthetic photometry as a proof 
of principle. We point out that MAP has the potential to constrain the eccentricity, 
and thus habitability, of Earth-like planets in the absence of radial velocity data, which 
is likely for terrestrial-mass objects. 

Key words: planets and satellites: general — eclipses — methods: numerical — 
planetary systems — techniques: photometric 



1 INTRODUCTION 

In February of 2011, 1235 Kepler transiting candidate plan- 
ets were announced by Borucki et al. (2011), amongst which 
the majority are expected to be genuine (Morton & John- 
son 2011). At the latest counting, the score has since risen to 
1781 (Rowe et al. 2011) and is expected to continue rising. 
Due to the unprecedented yield of new transiting planet can- 
didates, follow-up with radial velocity (RV) measurements 
is generally not feasible due to both the typical faintness 
of the targets and the intensive nature of the required tele- 
scope time for so many targets. Historically, radial velocity 
has emerged as the tool of choice to confirm transiting can- 



didates and so the Kepler team have devoted considerable 
effort to find ways to confirm candidates without the need 
for RV. This has led to some pioneering techniques such as 
blend analysis (Torres et al. 2011; Fressin et al. 2011) and 
confirmation through transit timing variations (TTV) (Hol- 
man et al. 2011; Lissauer et al. 2011). 

However, even though transiting candidates have been 
shown to be confirmable without RV, its absence means that 
the orbital eccentricity, e, of the planets cannot be deter- 
mined (unless very strong TTVs are detected). One remain- 
ing avenue to constrain e is to detect an occultation event. 
Occultations occur exactly half an orbital period after the 
transit event for a circular orbit ^ but become offset for eccen- 
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There also exists a small light travel time across the system 
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trie orbits, thus offering a potential diagnostic of eccentricity 
(details on the precise obtainable constraints are provided 
in Kipping 2011). Such occultations are due to a combinar- 
tion of reflected light from the planet and thermal emis- 
sion (which is very small in Kepler's visible bandpass e.g. 
Kipping & Spiegel 2011). Despite Kepler's ground-breaking 
photometric precision, most transiting planets discovered so 
far have not exhibited such events (one notable exception 
is the high albedo planet Kepler-7b, see Kipping & Bakos 
2011a). Therefore, in the majority of cases we are left with- 
out any way of characterizing the orbital eccentricity. 

Our principal motivation for addressing this problem 
stems from the fact that based upon the estimated fre- 
quency of Earth-like planets (Howard et al. 2010; Catan- 
zarite & Shao 2011; Wittenmyer et al. 2011), it seems prob- 
able that Kepler will detect numerous habitable-zone Earth- 
radius planets. After the initial detection, the natural ques- 
tion wc will be "can this world sustain life?" . A high orbital 
eccentricity causes a planet to spend the majority of its or- 
bit outside the habitable-zone and thus leads to potentially 
marginal or transient habitability (Williams & Pollard 2002; 
Dressing et al. 2010). For example, eccentricities in excess of 
0.3 for a habitable-zone planet around a Sun-like star cause 
equilibrium temperatures to vary by more than 100 K. The 
ability to measure, or at least constrain, eccentricity there- 
fore would greatly benefit the assessment of an exoplanet's 
habitability. 

Recently, Moorhcad et al. (2011) (Mil from here on in) 
have studied the problem in an effort to glean some infor- 
mation about the eccentricity of a transiting planet. Their 
approach is that if one knows the stellar density a-priori, say 
from stellar spectroscopy combined with evolution models, 
then one can compute the maximum allowed transit dura- 
tion for a planet on a circular orbit using (note that this 
approach was originally outlined in Ford et al. 2008): 



R^P 



(1) 



where T is the duration between the planet's centre 
crossing the stellar limb to exiting under the same condition, 
i?* is the stellar radius, P is the planet's orbital period and a 
is the planet's semi-major axis. The above equation assumes 
an equatorial transit and hence is the maximum duration 
possible. Mil discuss how if the observed transit duration 
exceeds this quantity (i.e. T/T^l^ > 1), this indicates that 
the orbit must be eccentric. The two weaknesses in this ap- 
proach are that: 1) The spcctroscopically determined value 
of R, has both a large statistical uncertainty and a large 
and unknown systematic uncertainty (for example. Brown 
et al. (2011) state that the Kepler Input Catalogue effective 
temperatures and radii estimations are reliable for Sun-like 
stars, but are "untrustworthy" for stars with T^s < 3750 K) . 
2) A planet can be eccentric yet still cause T /T^l^ < 1 i.e. 
only planets transiting near to apoapsc, the slowest part 
of the orbit, will be identified as eccentric, which is for 
w ~ 270° . The first weakness means the technique is model- 
dependent and that an individual system may not be re- 
liable due to possible systematic errors in Rt. The second 
weakness limits the scope of application of the technique to 
planets with 6 ~ and w ~ 270° , which it should be noted 



is the least probable value of oj from geometric priors (Kane 
& von Braun 2008). 

Mil show that the J?* uncertainty weakness can be 
overcome by adopting a statistical perspective. Even though 
an individual system may not be reliable, the bulk of systems 
should be and so any overall distributions which emerge 
should be reliable. The valuable technique of Mil allows 
us to actually say something about the eccentricities of the 
Kepler candidates as a whole. 

But what about individual systems? Or those which 
don't both fortuitously transit near (jj ~ 270° and have 
6 ~ 0? If the star has been studied with asteroseismology, 
then this prior can solve the riddle (as pointed out in Mil). 
This is scenario is discussed in more detail later in § 5.4. 
However, to date, relatively few targets have had asteroseis- 
mology studies conducted. We are therefore left with the 
quandary that it is usually impossible to say anything em- 
pirical about the orbital eccentricity of a transiting planet 
through photometry alone^. 



2 MAP: ANALYTIC CONSTRAINTS (a-MAP) 

2.1 Multiple Transiting Planet Systems 

We here present a solution for solving this riddle, which is 
applicable for multiple transiting planet systems. Although 
this may seem limited in application, in fact nearly 50% of 
all transiting planet candidates discovered by Kepler are in 
multi-planet systems (Rowc et al. 2011). Table 1 provides 
the most recent statistics on multiple-planet systems from 
Kepler. With this point established, we will now outline our 
proposed method. 

2.2 Light Curve Derived Stellar Density 

When a planet transits a star, consecutive transits provide 
the orbital period, P, and the light curve morphology con- 
tains information about the semi-major axis scaled in units 
of the stellar radius, a/R» (Seager & Mallen-Ornelas 2003; 
Kipping 2010). If it were possible to get a directly, then we 
use Kepler's Third Law to determine M* (we here assume 
Mp < M.): 



M. = 



GP2 



(2) 



Unfortunately, this is not possible. The light curve only 
lets us measure a/R*. Adding this into the above equation 
means we get the following: 



Ri ~ 
p* = 



GP-^Rl 
3ir{a/R*f 
GP2 



(3) 



Therefore, we can determine the stellar density from 
the transit light curve. This well-known trick, first pointed 
out by Seager & Mallen-Ornelas (2003), has been a powerful 
instrument in the toolbox of the exoplanetary scientist. It is 
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Table 1. Statistics regarding the number of planetary candidates found in multiple systems by Kepler. In this table, all examples of 
the word "planet" refer to a planetary candidate. Also, a "multi-system" refers to a solar system containing more than one transiting 
planet candidate. 



Si al isl ic 


.As iA Ijoriicki v\ al 


. {2011) 


.Vs oi iMiwe CI al. (zUiiJ 


Number or planets 






1 ^Ql 


Number of systems with planets 


997 




1296 


Number of multi-systems 


170 




328 


Number of single-systems 


827 




968 


% of systems which are multi-systems 


17.1% 




25.3% 


% of systems which axe single-systems 


82.9% 




74.7% 


Number of planets in a multi-system 


408 




813 


Number of planets in a single-system 


827 




968 


% of planets in a multi-system 


33.0% 




45.6% 


% of planets in a single-system 


67.0% 




54.4% 


2 planet systems 


115 




218 


3 planet systems 


45 




75 


4 planet systems 


8 




25 


5 planet systems 


1 




8 


6 ])laiiot svstoiiis 


1 




2 



common practice to determine p, from the light curve and 
then use stellar evolution models to estimate M* and 
separately. 

2.3 Eccentric Orbits 

As we discussed earher, without RV or an occultation we 
have no way of knowing what the orbital eccentricity of a 

transiting planet is. The major problem with this is that 
the determined value of a/f?, is heavily affected by orbital 
eccentricity. If we assume a circular orbit, but the orbit is 
really eccentric, then the determined value of a/R* would 
erroneously be (see Kipping 2010 for proof): 



(a/R,f 



1 -|- e sin w 



(a/R^) 



(4) 



In other words, 
(1 + esinw)/-\/]~ 



a/i?* is wrong by a factor given by 
e^. Note that the above equation is an 
approximate formula based upon the T°°* expression for 
the transit duration (Kipping 2010). The reliability of this 
approximate expression will be dealt with later in §2.7. 

If we determine a biased value for a/R*, then one will 
determine a biased value for p, too, since p* depends upon 
a/Rt, as seen in Equation 3. This means that the derived 
stellar density becomes 



where 



(1 + esincj)^ 
(l-e2)3/2 



(5) 



(6) 



and hence no reason to assume otherwise. For each planet, 
one may empirically determine pl^^k'- 

P?,? = *2P* (7) 

Since the parent stax is the same for planets, then one 
can divide these two equations to give: 



circ 



*1 
*2 



1 + ei sinwi 

1-1-62 sin IjJ2 



3/2 



(8) 



^1 , 

Note, that a similar equation to this appears in 
Ragozzine & Holman (2010), Equation 3, although the equsr 
tion is known to contain an error (D. Ragozzine personal 
communication). Nevertheless, the authors hint at the pos- 
sibility that it may possible to determine eccentricities via 
this ratio. 

We stress here that the term on the left-hand-side of 
Equation 8 is an observable. One can therefore see that it 

is possible to glean some information about the eccentric- 
ity of the system. However, the current form is not very 
informative. We have one observable, (p*'T/p?,2")i and four 
unknowns: ei, sinwi, 62 and sinci;2. 

Ultimately, the quantity which will affect potential hab- 
itability of a planet will be ek and not Wk- Therefore, the uik 
quantity is of lower interest to us. The sin w/c terms must lie 
in the range — 1 ^ sinwfc ^ which means we can con- 
struct a lower and upper bound inequality for the quantity 
612: 



2.4 Double Transiting Systems 

Consider two planets, dubbed with subscripts "1" and "2", 
which have been observed to transit the same star. Let us 
assume that we fit these light curves assuming ei = and 
62 = 0, since we have no RVs, occultation or strong TTVs 



Gi 



Oi 



„circ 

P*,i 



2/3 



1 -I- ei sinoji 

1 + 62 sin UJ2 



(9) 
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Which must satisfy: 

(lTt)(^)<«'-(S)(^) 

One may now replace €2 = €2161 where £21 = 62/61 and 
expand to first-order in ei: 

-(ei + 62) + 0[ei]^ < ^ (ei + e2) + 0[eif 

^ 012 — 1 ,^ 

.•.6l+62> (11) 

If one instead replaces 61 = €1262 in Equation 10 and 
expands to first order in 62, an identical equation is obtained. 

Before moving onto triple systems, we discuss a final 
subtlety. The fraction Q12 can be inverted to give O21 and 
yet the same inequality is derivable as Equation 10, namely 
we have: 

Following the same steps as before, means we arrive at: 

, 621 — 1 , , 

61 + 62 > (13) 

So wo have two equations, both of which must hold, 
for constraining the (ei + 62) combination. Both Equa- 
tion 11&13 require that Qij > 1 in order to place positive 
constraints on this value (a negative value has no physical 
meaning). Clearly, if O12 < 1 then 02i > 1 and vice versa. 
Thus we can always construct a physically meaningful con- 
straint on (ei + 62) by taking the maximum of the two. 

In practice, we wish to produce a posterior distribution 
of (ei + 62) based upon the posterior of O12 or O21. We 
can choose to use either version of Qij but not both i.e. we 
cannot create a posterior which swaps between the two Qij 
versions. The simplest thing is to produce two posteriors 
and then select the one which provides the most meaning- 
ful constraints. This selection can be done visually, or by 
say taking the median of both posteriors and choosing the 
largest. 

2.5 Reliability of the a-MAP Inequality 

Equations 11&13 are approximate equations valid to first 
order in Ck only. Therefore, the reliability of the inequality 
will deteriorate for large efc. To test the accuracy of these 
expressions, we generated some random values for ei, uji, 32 
and 0J2- The oJk values have uniform distributions between 
and 27r and the Ck values have uniform distributions between 
and 6max. We generated these random values 10® times and 
tested if the inequality in Equation 11 was true or not each 
time (the accuracy of Equation 11 will be the same as that of 
Equation 13 due to symmetry arguments). As an example, 
using 6max = 0.25, the inequality is true in 91.9% of all 
of the Monte Carlo simulations. In Figure 1, we show the 
percentage of trials for which the inequality is correct as a 
function of emax, which reveals that the inequality provides 
useful eccentricity constraints in the absence of any other 
information and is ^90% reliable for 6max ^ 0.30. 




0.0 0.2 0.4 0.6 0.8 1.0 



Maximum eccentricity (of either planet) 

Figure 1. The 'rcliabH/itij of llie. i/ii equality in Equation 11 as a 
function of the maxiniuni allowed eceentricity, assuming a uni- 
form distribution in ej. (dashed line) and also a non-uniform 
physically motivated distribution (solid line). The reliability is 
^90% for efc ^ 0.30 for the uniform case, and e^ ^ 0.65 for the 
non-uniform case. 



We also tried using a potentially more realistic non- 
uniform distribution eccentricity distribution using a mix- 
ture of an exponential and a Rayleigh distribution (see Juric 
& Tremaine 2008; Zakamska et al. 2011): 

P(6fc) = aAexp[-A6fc] + {1 - q)^ exp[-e^/(2<Te'] (14) 

The values of the constants were found by fitting the 
distribution of eccentricities in known multi-planet systems 
measured from radial velocity surveys using only systems 
with measured eccentricities, which find a = 0.38, A = 15 
and fTe = 0.17 (Steffen et al. 2010). Finally, this distribu- 
tion can produce values of greater than unity, and so we 
ignored any simulations where ep > emax for cither planet. 
Using emax = 1, we found that 87.0% of simulations agreed 
with the inequality presented in Equation 11, and ^90% 
agree for emax ^ 0.65 (Figure 1 shows dependency of this 
percentage with emax). 



2.6 Triple Transiting Systems 

For three-planet systems, one may construct three Qij ratios 
meaning we now have six unknowns and three observables. 

By blanketing out the ujk terms in a similar way described 
for double-transiting systems (and thus limiting ourselves to 
lower bounds on e^) we are left with three unknowns (ei, 62 
& 63) and three observables (612, 623 & G31, following a nu- 
merical cyclic). One would therefore presume that it should 
be possible to constrain the eccentricities individually rather 
than limiting ourselves to combination terms. In total, we 
have three inequalities, in analogy to the double-transiting 
case: 



Mulitbody Asterodensity Profiling (MAP) 5 



ei + 62 > -{612 -1,621-1} 

62 + 63 > ^{623 - 1,632 - 1} 

63 + 61 >i{631- 1,613-1} 



(15) 

(16) 
(17) 



One may naively assume that these may converted into 
individual limits by solving simultaneously. However, such 
an operation requires subtracting one inequality from an- 
other, which is a strictly forbidden operation. With this bar, 
it is not possible to solve the expressions and so the furthest 
we can ever take a-MAP is to produce the inequalities of 
Equations 11&13. Naturally, this leads us to consider alter- 
native methods which are not based upon analytic methods. 
However, before we do, we will pause to evaluate the range 
of parameters under which Equation 5 (the ^'-equation) is 
valid. This is crucial since all of the a-MAP expressions are 
built around this equation and even the later numerical tech- 
niques have the same dependency. 



2.7 Reliable Parameter Range for the ^-Equation 

Equation 5 has been revealed to be the key to unlocking 
some information about the eccentricity of transiting plan- 
ets. It is worth, though, pausing to evaluate the reliability of 
this expression, since the equation is an approximation, as 
explicitly stated in Kipping (2010). It should also be stressed 
that the non-analytic method discussed later (§3) relies on 
the simplicity of the ^'-equation too, and would not work 
using a more elaborate form. The reliable range of the 'I' 
equation informs the reliable range of our proposed tech- 
nique in general, and thus is crucial for understanding the 
range of applicability of our new method. 

The Vt-equation describes how erroneous the derived 
stellar density would be if we assumed a circular orbit for 
an eccentric planet. As was seen in Equation 3, the erro- 
neous pt value really comes from an erroneous a/Rt value. 
Since the derived orbital period will be reliable irrespective 
of the orbital eccentricity, the only cause of p, deviating 
from the true value is because (a/7?«)^ deviates from the 
true value. Using a set of approximate expressions for the 
transit duration, which are accurate to 99.9% for planets of 
|esinaj| < 0.5 and \ecosio\ < 0.85, Kipping (2010) showed 
that if one assumes a circular orbit then the derived value 
(a/Ri.) differs from the true value via: 



[{a/R.y"^f = ^V+g'^^^/oT + I^""]' (18) 



sin2(Ti,47r/P) 



2 



sm 



: arcsin( 



-|- sin 



„2 



: arcsin( 



^(l_p)2_b2 

{a/R*)Qcami 

^(l+p)2_62. 



{a/R.,)^ 



{a/R,)Qc sini 



: arcsin( 



^(l+p)2_fc2 



(19) 



where Qc = (1 — e^) / [1 + esinu) and p is the ratio-of- 
radii. The relative difference between our approximate ex- 
pression for the stellar density (i.e. the ^i-equation) and the 
more accurate value is therefore given by the LHS of the 
following: 



ia/R,f^ 



[{a/R.Y 



< t 



(20) 



Where t is the tolerance level for the desired level of 
accuracy. For example, a typical choice might he t = 10~^ 
indicating 99.9% accuracy in the ^-equation. For brevity, 
we do not write out the full form of the above expression. It 
is trivial to show that it is maximized for w = 7r/2 and 6 = 1 
and therefore if we satisfy a given tolerance level under these 
conditions then we can be sure the equation is always valid. 
Eliminating these two terms accordingly, one may then take 
the limit of the resulting equation for when p — ^ 0, corre- 
sponding to the small-planet approximation which is essen- 
tially valid for p < 0.1, encompassing almost all transiting 
planets. This leads to the far simpler expression: 



(a/R^' 



(1 



(1 - e2)3/2 



{a/R*f{l - efjl + e) - e(4 - 3e -h 
(l-e)3 



-3/2 



>l-t 

(21) 



where 6'^"^'^ is given by: 



This expression indicates that our accuracy becomes 
worst for low (a/R*) and high e values. Numerically solving 
for the maximum allowed e as a function of (a/R*) may be 
accomplished for a given t level to illustrate some typical 
constraints. In Figure 2, we show the case for t = 10~^ 
(dashed) and t = 10"'^ (dotted). 

Since even the highest precision measurement uncer- 
tainty on p« is around 1% (e.g. Kipping & Bakos 2011b), 
t = 10~^ is sufficient for our purposes. For some typical 
values of (a/7?,) of (a/i?.) =10, 100, 300 & 1000 (corre- 
sponding to a hot-, warm-, temperate- and cold-planet re- 
spectively, assuming a Solar-star) we find Cmax =0.15, 0.88, 
0.96 & 0.99 respectively. If these limits are exceeded, the 
techniques presented later in this paper will still detect or- 
bital eccentricity but the actual determination of the e will 
obviously be subject to a systematic error. 

This limitation may seem to be a significant disadvan- 
tage but in reality one does not expect to find multiple 
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1.0 




a/«. 



Figure 2. Maximum allowed orbital eccentricity versus (a/Rt) 
in order for the -equation to remain a valid approximation. The 
two lines correspond to two accuracy tolerance levels. The dashed 
line is for 99% accuracy and the solid is for 99.9%. Small humps 
are due to numerical errors. 



planet systems with very high eccentricities. The most fear 
sible scenario where the expressions would become invalid 
is for a multiple planet system featuring one very close-in 
planet. However, in such a case, one would not expect the 
planet to be highly eccentric anyway due to rapid tidal cir- 
cularization at such distances. Nevertheless, if the planet is 
both close-in and highly eccentric, one can simply discount 
the planet during the MAP analysis. 



3 MAP: NUMERICAL CONSTRAINTS 
(n-MAP) 

3.1 Description of the Algorithm 

Whilst the approximate inequalities derived in §2 are useful, 
they do not take advantage of the extra information offered 
by > 2 planets and are invalid at large eccentricities since 
we are ignoring terms of 0[e\] or larger. To address this, 
one may consider a more brute-force approach through nu- 
merical methods. 

§2 has shown that converting Qij — > {e/bjWfc} is degen- 
erate and in general we can only obtain some lower bounds 
on Cfe. However, it is trivial to perform the reverse operation 
and convert {ek,u)k} -> ©ij- Therefore, although compu- 
tationally demanding, one could create a multi-dimensional 
grid of all possible and ujk values, for all planets, and 
compute the Qij values at each grid point. The Qij value at 
each grid point could then be compared to the observed val- 
ues of Qij to evaluate the likelihood of the {e/c,aj/b} vector at 
the specified grid location. This likelihood can then be used 
to create maps of the permitted/excluded parameter regions 
for {efe,cjfc}. This numerical approach is an extension to the 
analytic approximations used earlier. Both techniques can 
be regarded as what we call "Multibody Asterodensity Pro- 
filing" (MAP). We distinguish between the two approaches 
as analytic-MAP (a-MAP) and numerical- MAP (n-MAP). 

Although we have just described a grid search in 
the previous paragraph, a much more computationally 
efficient and powerful numerical technique is to adopt 



a Markov Chain Monte Carlo (MCMC) algorithm. Our 
overall approach to inunerical-MAP utilizes two-stages of 
a Metropolis-Hastings MCMC fitting routine. The first is 
for the light curve fits and the second is performed once 
the light curve MCMCs are complete. Experimentation 
with combining the two stages into one larger MCMC 
yielded inordinately large computation times, whereas the 
two-stage technique provides results in just a few hours on 
a typical workstation. Our method can be loosely described 
via the following algorithm: 

Stage 1 

1.1 For an n-body system, fit transit light curves for each 
planet independently, assuming a circular orbit, using a 
Markov Chain Monte Carlo (MCMC) algorithm with the 
Metropolis-Hastings rule. 

1.2 Compute n!/[2!(n — 2)!] marginalized posterior distri- 
butions for Qij for i = 1,2, ...n, j = 1, 2, ...n and i ^ j. Each 
posterior is ensured to have A down-sampled points from a 
well-mixed and converged chain. 

1.3 Normalize the posteriors such that they become equal 
to the probability density function (PDF) of each Qij . We 
dub these PDFs as PBF{Qij). 

Stage 2 

2.1 Define a starting point for a new MCMC chain with 
a fitting parameter set A;, = {efc,6,a;fc,6} for k = 1,2, ...n, 
where b is understood to represent the b*^ accepted MCMC 
trial. 

2.2 Using Equation 9, evaluate Qij,b for i = 1,2, ...n, j = 
1, 2, ...n and i ^ j. 

2.3a Ifthe trial value of Oij^fc > Mode[PDF(0ij)], then count 
the number of realizations in the PDF(Oij) which fall in the 
range of Mode[PDF(0ij)] < Qij < Qij^b and define this 
integer as rriij^b- 

2.3b If the trial value of Qij,b < Mode[PDF(eij )], then count 
the number of realizations in the PDF(Gij) which fall in the 
range of Qij^b < Qij < Mode[PDF(Oij)] and define this 
integer as rriij^b. 

2.4 Define the of the 6* MCMC trial as Y^i^j Xij,b where 

X?,,6 = {V2Err'[2rn,j,b/A]f. 

2.5 Accept/Reject trial point following the Metropolis- 
Hastings rule and loop the MCMC in the usual manner until 
B trials have been accepted. 

By the end of the algorithm, we have obtained a joint- 
posterior for the {ek,iOk} vector revealing those regions of 
parameter space which are excluded and those which are 
more probable. The merit function (inverse of the likelihood) 
is thus given by: 

2 \ ^ / Qij,obB ~ ©ij, model | /OON 
^--=1^1^ ^) j ^''^ 

where it is understood that a is determined by numer- 
ically integrating the probability density function of ©ij,obs 
(steps 2.3&2.4), rather than simply counting the number 
of la error bars between the model and observed value of 
Qij. The advantage of doing this is that we are able to fully 
account for non-Gaussian ©ij,obs posteriors, which are com- 
mon as will become evident in §4.1, 4.2 & A3. The disadvan- 
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tage is significantly increased computation time since each 
MCMC trial of n-MAP requires ~ n!/[2!(n - 2)!]yl evalua- 
tions. 

We also point out that one may add additional infor- 
mation into the technique at this stage. For example, if ra- 
dial velocity or TTV information exists and can be used to 
place further constraints on Ck, then additional components 
to the total merit function cam suitably appended. However, 
for the remainder of this work we focus on MAP alone to 
demonstrate the use of this technique as a unique type of 
observable. 

To compute the inverse error function, we use the ap- 
proximation of Winitzk (2006), which is accurate to 4x 10~^ 
over the interval < a; < 1: 



Erf"^(x-) 



2 



log(l - x^) 



2_ 



log(l 



^log(l 



X2) 



1/2 



(23) 



where P = [8(7r — 3)]/[37r(4 — tt)]. In the next subsection 
we discuss how we enforce a uniform prior in (p«'fc^)^^^, which 
enables a uniform prior in the Qij posteriors. 



3.2 Light Curve Fitting Parameters 

The choice of fitting parameters for the transit light curve 
has a broad diversity within the exoplanet literature and 
yet a significant impact on the derived results and efficacy 
of a light curve fitting algorithm. In this sxibsection, wc will 
describe the light curve fitting parameter set which yields 
the most reliable results for the specific purposes of MAP. 
In order to accomplish this, we will briefly overview the basic 
properties of the light curve. 



3.2.1 Understanding the trapezoid light curve 

The transit light curve is well-approximated by a trapezoid 
in the limit of negligible limb darkening. A trapezoid is de- 
scribed by four parameters only: T14, the duration from the 
l''*-to-4* contact; T23, the duration from the 2"'^-to-3''* con- 
tact; 5, the depth of the trapezoid and tmid, the mid-point 
of the trapezoid in time. The parameter set may be written 
as r = {5, ri4, r23, tmid}. It is easy to see that one could 
alternatively use the ingress/egress duration {T12) and the 
full- width-half-maximum duration (W) instead of T14 and 
T23. 



3.2.2 Understanding the circular orbit light curve 

Consider a transiting planet on a circular orbit but with 
limb darkening on the star. It is easy to show that four pa- 
rameters only can still be used to completely describe the 
light curve, even though the light curve is no longer mor- 
phologically trapezoidal (Kipping 2011). These terms are: 
(a/7?,), the semi-major axis of the planetary orbit around 
the star in units of the stellar radius; 6, the sky-projected 
distance between the planet and the star at the instant 
of inferior conjunction in units of the stellar radius (often 



called the impact parameter); p^, the square of the ratio of 
the planet's radius to the stellar radius and r, the instant 
when the sky-projected planet-star separation is minimized 
in proximity to the instant of inferior conjunction. So we 
have f = {a/R,),b,T}. 

One can easily appreciate that p^ replaces 5 (but are 
equivalent for a non-limb darkened star) and r replaces imid 
(but are equivalent for circular orbits). Further, {{a/R»),b} 
replace {714,123}. However,as shown by Seager & Mallen- 
Ornelas (2003), these terms are interchangeable via: 



lim b ■■ 

e->0 



lim(a/Ji*) = 



2 Bin^(T23 7r/P) 
sin^lTllTT/P) 



sin- (1-2 



a^(Ti4,r/P) 
{l+pf-[b''{l-Sm\Tl47T/P) 



1/2 



1/2 



sin2(r237r/P) 



(24) 



(25) 



Therefore, one has the choice as to whether one uses 
{{a/Rt),b} or {ri4,r23}. Indeed, one can also legitimately 
use many other combinations which are interchangeable, 
such as {W,Ti2} (Carter et al. 2008), {f,Ti2} (Kipping 
2010), {{C/R*),b'^} (Bakos et al. 2007), {{r/R,),b'^} (Kip- 
ping 2010), etc. This already raises the question as to what 
parameter set should be used. The two terms are problem- 
atic in that they typically exhibit mutual correlation and so 
care must be taken in their selection. 



3.2.3 Understanding the eccentric orbit light curve 

For an eccentric orbit, the morphology of the light curve is 
essentially unchanged. The signal of asymmetry is negligible 
and will rarely affect measurements for even extreme cases 
(Kipping 2008; Winn 2010). As a result, the eccentric terms 
e and uj (eccentricity and position of pericentre) are hid- 
den from view and cannot be determined by simply fitting 
a light curve (obviously, for multi-planet systems an alter- 
native, more subtle strategy exists in the form of MAP). 
Consequently, the same parameter set applies for eccentric 
orbits as for circular orbits i.e. one can use, for example, 
r = {p^,{a/R*),b,T}. The only difference is that the we 
have to declare values for e and u} during the fits. These ec- 
centric terms do affect the relationship between {{a/R,),b} 
and duration related terms and a modified form of Equa- 
tion 25 should be used, as presented in Kipping (2010). In- 
deed, it is these differences which fundamentally allow MAP 
to work. 



3.2.4 Choosing a parameter set 

In this work, the term which we are interested in is the de- 
rived value of Qij, which has been established to contain 

information about the orbital eccentricity. This term is sim- 
ply the ratio of {pl"'^)^'''^ values. As is discussed in §3.4, 
uniform priors in the eccentricity terms for n-MAP can be 
implemented by ensuring a uniform prior in (p"^'^)'^^^. By 
fitting for, say {p^ , (a/Rt,), 6, r^}, we necessarily assume uni- 
form priors on those terms. However, since {pl"'^)'^^^ is an 
intricate function of these terms, it will not have a uniform 
prior. 

A simple but effective solution to this problem is to fit 
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for {pT")'^^'-^ directly. This term may be easily converted to 
(a/Ri.) via: 



{a/R* 



GP^pT 



1/3 



(26) 



This leaves one of the two problematic terms assigned, 
but still leaves us with options for the other. For example, 
should we use {p' ,{pt'^f/\b,n} , or {p^ ^12, r^}, 

or {p^, (pt"^'^)^''^, cos(i), Ti} or many other possible permutsr 
tions? 



3.2.5 The other term 

Choosing this term is non-trivial. If for example, we chose 
6, we would be faced with the issue that 6 < is unphysical 
and thus a boundary condition exists at 6 = 0. In a Marko- 
vian sense, jumps to negative 6 values are rejected thus re- 
sulting in a disequilibrium between the number of positive 
and negative jumps. This in turn means that the h posterior 
will be biased and overestimate the true value. Since inter- 
parameter correlations exist between the light curve fitting 
parameters, a bias in b induces a bias in (p'^^^^Y^''' . Terms 
being correlated in itself is not a major problem, it slows 
down our algorithm but an accurate result can still be ob- 
tained with a sufficient number of trials. However, if one of 
the correlated terms is biased then the fact the terms arc 
correlated to one another becomes a problem, since now all 
terms will become biased. 

The boundary condition in h in therefore a serious is- 
sue. A similar situation is well-known to exist for e with 
radial velocity fits (Lucy &: Sweeney 1971). One could pro- 
pose that using a duration-based term such as ri4, T or T\2 
would avert such a problem. However, for certain duration 
jumps, the derived impact parameter still falls out as being 
unphysical (in this case imaginary). These unphysical trials 
can bo discarded but that again introduces a bias. There- 
fore, any other duration related parameter would also be a 
boundary-condition-limited parameter. 

To allow h to maJfe Markovian steps, we let h go to 
negative values. Since the light curve and duration-related 
terms are always generated using 6^, then the physicality 
of 6 < solution is irrelevant mathematically speaking. We 
found that this yielded solutions consistent with test cases 
and thus seems to solve the problem. 

In addition to {p^ ,{pT''f''\b,n}, we also allow each 
transit epoch to have a unique out-of-transit normalization 
factor, OOTjn where m denotes the epoch number. The ze- 
roth epoch is defined to be that which has the lowest mutual 
correlation to the orbital period. Finally, the orbital period 
is a free parameter too. 



3.3 Direct n-MAP Priors 

In most cases (system with less than 
ets), there are more free parameters in 
servables and so the problem is under-i 
unique solution. Despite this, contours 
may still be computed through Monte 
our case, the Monte Carlo technique 
Chain Monte Carlo (MCMC) with the 



five transiting plan- 
the model than ob- 
constrained with no 
of the error surface 
Carlo techniques. In 
of choice is Markov 
Metropolis-Hastings 



algorithm. We note that an analogous situation arises for in- 
terpreting the atmospheres of exoplanets where atmosphere 
models tend to have more free parameters than the num- 
ber of observations. In this example, a similar solution as 
n-MAP has been adopted in works such as Madhusudhan 
& Scagcr (2009) and Madhusudhan et al. (2011). A direct 
comparison of n-MAP to the methods of Madhusudhan et 
al. (2011) is discussed in §3.5. 

We stress here that because the problem is under- 
constrained (except for n ^ 5), the results will be strongly 
affected by the choice of priors. This is in contrast to a highly 
constrained problem where the data drive the result to the 
same solution with only a minor dependency upon the choice 
of priors. Thus, the choice of priors can be understood to 
have a significant impact on the derived results (see Ford et 
al. 2005 for a detailed discussion on the effect of priors when 
fitting exoplanet data). The results should consequently be 
always quoted in unison with the adopted prior used to infer 
them. 

These ideas are more formally expressed through Bayes' 
theorem. Let us denote the eccentricity parameters which we 
fit for in n-MAP as A = {efc, Wfc, ...} for k = l,n. We further 
use M. to represent the model and V to represent the data: 



V(K\V,M) = 



P(A|.M)P(P|A,.M) 
fP{A\M)P{V\A,M)dA 



(27) 



In our case, the "data" is the observed ratios of 
(pi^i/pfjf^'^, denoted by the term 6,^: 



P(A|e,M) 



P(A|.Vl)P(e|A,.Vl) 
/P(A|A4)P(e|A,A^)dA 
P(A|.A/l)P(e|A,>f) 

p(e) 



(28) 



Given that the problem is under-constrained, clearly 
the choice of this prior will have a significant impact on the 
derived joint probability distribution from n-MAP. There 
are two plausible paths to adopt: 

■ Assume complete ignorance for the a-priori knowledge 

of A 

■ Adopt a prior based upon dynamics and/or known prior 
distribution of eccentricity 

For the former, complete ignorance can be easily im- 
plemented by adopting a uniform prior in A. This would 
take the form of a uniform prior between < ej, < 1 and 
^ a;fe < 2tt. One advantage of this choice is that any re- 
sults derived from n-MAP can be understood to be directly 
due to the MAP technique rather than any prior biases. 
One disadvantage is that we know that a system of multiple 
planets is unlikely to survive with high eccentricities and we 
arc essentially ignoring this fact. However, this could also be 
considered a potential advantage in that a system where n- 
MAP strongly prefers a dynamically unstable solution may 
indicate that the system is in fact a false positive. 

For the second option, a typical procedure is outlined 
in the recent work of Stcffcn ct al. (2010) for five candidate 
multiple transiting planet systems detected by Kepler. Here, 
the authors adopted a prior distribution in based upon 
the same distribution discussed earlier in §2.5 i.e. a mixture 
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of an exponential and a Rayleigh distribution following Juric 
& Tremaine (2008) and Zakamska et al. (2011) and provided 
earlier in Equation 14. 

As already touched on, in this work we prefer to present 
n-MAP results with uniform priors in A (i.e. assume total 
ignorance) for the sake of demonstrating this new technique, 
but future works could make use of more sophisticated priors 
like those of StefFen et al. (2010). By making this choice, the 
derived constraints in this work can be understood to be 
completely due to the n-MAP method alone, and not due to 
the impact of priors. 

With this choice, it is arguably better to think about the 
n-MAP results in terms of "allowability-space" rather than 
probability space. An area of high-density from n-MAP sig- 
nifies a region where lots of combinations of parameters can 
reproduce pl^^k'^ which agree well with the data. An area of 
low-density signifies a region where very few combinations of 
parameters can reproduce s which agree with the data. 
An area of null-density signifies that absolutely no combi- 
nation of parameters can reproduce the observed p^^^'s. 



3.4 Indirect n-MAP Priors 

Aside from the priors in the MCMC chain of the n-MAP 
phase, priors also affect n-MAP indirectly via the light curve 
fits. If we, for example, fit for pf^ in the light curve fits, the 
prior on Qij will be non-uniform. Even if 0jj has a uniform 
prior, it is not immediately obvious that this will translate 
to uniform priors in and ujk- 

Consider first that one executes n-MAP with the Oij 
terms behaving as uniformly distributed parameters. Qij is 
used to compute the of the n-MAP realizations in step 
2.3 (as described in §3.1). This calculation requires that we 
know the mode of Oij, which is a meaningless concept for a 
uniform distribution. Nevertheless, one can easily appreciate 
that the must be the same for all MCMC realizations 
of {ek,ujk}, since Qij is uniformly distributed. In therefore 
follows that all MCMC realizations will be accepted under 
the Metropolis-Hastings rule. If all trials arc accepted, then 
this is equivalent to the case KOI-SOP described in §A1, 
which simply reproduces the behaviour of the direct priors 
i.e. uniform in {ek,u>k}- Consequently, uniform priors in Qij 
are something to be desired since it does not cause any bias 
in the resulting n-MAP procedure. 

With this point established, the next question is how 
can we ensure uniform priors are produced for Qij? Since 
Qij is the ratio of {pfA )^^^/{pf^^f^^, one first step would 
be to adopt uniform priors in {pl^^k)'^^'^ ■ However, the ratio 
of two uniform priors does not produce a uniform prior itself, 
rather we have: 



PiQij) = { 



1 

2 

1 

2e? 







< Gij < 1 

O^J > 1 

otherwise 



(29) 



Ergo, by fitting the transit light curves for {pl"k)^^^, 
we adopt a uniform prior in this term but the derived Qij 
terms will be non- uniform for Qij > 1. Our solution for 
tackling this is to adapt the n-MAP algorithm. The solution 
lies in the fact that Qij is uniform for values < 1. If we 
generate a {ek,ujk} realization which causes Qij > 1, we 



may simply use Qji instead, which must lie in the range 
< Qji < 1 and therefore must be uniformly distributed 
according to Equation 29. By implementing this condition, 
we ensure that both the direct and indirect priors in n-MAP 
are uniform in {ek,ujk}- 



3.5 Difference Between n-MAP and the 
Madhusudhan et al. (2011) Technique 

There are several methodological similarities between the 
n-MAP technique to constrain orbital eccentricity and the 
rmmcrical mctliods proposed by Madhusudhan & Seager 
(2009), and subsequent papers, to determine exoplanet at- 
mospheric composition. In these papers, the authors pro- 
duced parameter space maps showing the points which agree 
to the data to within x^ < 1, X^ < 2 and x^ < 3 to denote 
different error surfaces. As of Madhusudhan et al. (2011), 
MCMC methods were used for the parameter space explo- 
ration rather than grid methods which were used in Mad- 
husudhan & Seager (2009). Despite switching to MCMC, 
the presentation of results remained largely unchanged with 
X^ < 1, etc surfaces still being plotted. Such a presentation 
does not take advantage of the fact the MCMC technique 
inherently computes the probability density over the pa- 
rameter space, rather than merely outputting the likelihood 
of individual realizations. To accomplish this, one simply 
computes the parameter space regions in which the MCMC 
spends the majority of its time. These regions represent the 
high probability density areas. Therefore, we can see that 
there are two possible paths by which to proceed. 

The analogy between n-MAP and the method of Mad- 
husudhan et al. (2011) breaks down here. For us, the first 
way of presenting the results would be a far less useful di- 
agnostic of the parameter space. This is because almost any 
Bi-Bj combination can be found to give a x^ < 1 for the cor- 
rect tuning of uJi and u)j . However, the fine tuning of these 
parameters must be so precise, that very few MCMC real- 
izations find such values. Nevertheless, we typically compute 
O[10®] points in the MCMC chain and so these improba- 
ble locations will eventually be visited by the chain. Conse- 
quently, if we plotted the minimum x^ in a rasterized grid of 
ei-ej, analogous to the presentation in Madhusudhan et al. 
(2011), we would essentially find a region where any solution 
is permitted. In contrast, plotting the probability density re- 
gions in Bi-Bj space automatically accounts for the fact that 
despite these locations yielding a low x'^, they require very 
precise tunings of uji and ujj. This is the key difference in 
the presentation of our results. 

A way to visualize this more easily to consider the result 
we would obtain for a two planets on a circular orbit. An 
numerical example will illuminate this issue more fully. Let 
us say p^'.T = 1.00 ± 0.0071 gcm"'^ and pl%^ = 1.00 ±0.0071 
,gcm^'\ which is consistent with that which would be ob- 
tained for two planets on circular orbits. These derived stel- 
lar densities would yield 612 = 1.000±0.0067. Now consider 
that during the MCMC parameter exploration, one realiza- 
tion is attempted where ei = 0.9 and B2 = 0.4. This would 
seem to be a position that one would expect to be highly 
excluded by the n-MAP technique. We have: 
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Figure 3. For e\ = 0.9 and 62 = 0.4, yet ©12 = 1.00, the above 
shows the range of the imaginary component of 0J2 values which 
allow a solution of = in solid and = 9 in dotted, as a 
function of loi . Solutions requiring non-zero imaginary L02 can 
never be reached by the n-MAP algorithm. However, a narrow 
range of uii allows both real uji and UI2 solutions. For this rea- 
son, although improbable, n-MAP can find a low \ solution at 
virtually all Ci — ej coordinates. However, these states are visited 
with a low frequency due to the above constraints. 



ei2 = 



1 + ei sina;i 



1 + 62 sin W2 / \ 1 — ef 
Solving for G12 = 1 for U12, we have: 



(30) 



UI2 = arcsm 



-ei(l-ef) 



(31) 

Plotting the imaginary component of this equation as a 
function of wi for the fixed values of ei = 0.9 and 62 = 0.4 
(Figure 3), a narrow range of wi values can still yield a 
real result. Specifically, in this example, 17.2% of the al- 
lowed oji range can yield a solution. One solution occurs at 
wi = 217°. We can therefore imagine the MCMC algorithm 
landing upon this value too, thus permitting a physical so- 
lution. Even if it happens to land within this 17.2% sliver 
of LOi, we also require a very narrow range of UJ2. Here, we 
require uJi = (354.8 ± 0.5)° in order to obtain < 1 in n- 
MAP. This is a range of just 0.3% of all possible L02 values. 

Putting this all together, for a genuinely circular orbit 
system, a random realization of ei — 0.9 and 62 — 0.4 can 
still have some n-MAP trials of < 1- However, the frac- 
tion of times which n-MAP will succeed in doing this will be 
less than 0.05% (for uniform priors in u>k). Nevertheless, for 
long chain lengths we should expect n-MAP to find x^ < 1 at 
virtually every single point in Ci-Cj parameter space. There- 
fore, presenting plots of points where realizations yielded 
various x^ is of limited value. What is much more useful, 
is to plot the density of realizations across the Ci-ej param- 
eter space resulting from n-MAP. This inherently accounts 
for the fact certain Ci-ej combinations require fine tunings 
of the other terms in order to produce an acceptable trial. 
This is the chief difference between n-MAP and the method 
of Madhusudhan et al. (2011). 



3.6 MCMC Diagnostics 

We will later show applications of our algorithm to synthetic 
data sets. In these cases, it is important to ensure that the 
MCMC fits for both the transit light curves and n-MAP 
achieve i) adequate mixing ii) adequate convergence. 



3. 6. 1 Bum-in 

Before either of these diagnostics can be computed, it is im- 
portant to remove the pre-burn trials of the MCMC. These 
trials are highly dependent upon the initial starting point 
of the chain and thus it is important to burn-out the initial 
part of the chain. We will use the same strategy as Tegmark 
et al. (2004) for this. Specifically, we compute the median 
of all accepted MCMC trials and then burn-out the initial 
trials up to the point when the drops below the median 
value. Burn-out is typically very rapid and occurs within a 
few dozen trials. 



3.6.2 Mixing 

To determine whether our chains are sxifliciently mixed, we 
compute the effective length of the chain (see Tegmark et 
al. 2004). Each free parameter has its own unique effective 
length and so we always conservatively adopt the lowest 
effective length in reporting the final diagnostics. Broadly 
speaking, we wish to reach a point where the lowest effective 
length ^ 1, such that meaningful statistics can be inferred. 
In this work, we set the goal that the lowest effective length 
> 1000. This is sometimes achieved by combining multiple 
chains rather than simply extending the length of a single 
chain. The advantage of this is that the chains may be run 
simultaneously with parallel processing and yet still com- 
bined at the end provided that the burn-in trials have been 
removed and that each chain has reached adequate conver- 
gence. 



3. 6. 3 Convergence 

Convergence for each free parameter in each chain may be 
checked by computing the Geweke (1992) statistic. This sim- 
ple statistic compares a given parameter's value at the be- 
ginning of the chain and at the end of the chain, accounting 
for the variation due to parameter exploration. It is essen- 
tially characterizes the number of sigmas difference these 
two points. In a converged chain, we require that the Geweke 
(1992) statistic < 1 and certainly ^ 3. 

Convergence is not generally expected with n-MAP 
since the problem is usually under-constrained. Thus the 
application of the Geweke (1992) statistic is not employed 
for n-MAP results. 



3. 6. 4 Down- Sampling 

Sometimes a fit requires either a very long chain or mul- 
tiple chains which are combined. This can lead to a very 
large number of points in the final combined chain, of order 
10''-10^. In general, this many points is excessive to build 
reliable posteriors. One significant disadvantage of this is 
that n-MAP must count the number of trials below/above 
various Qij thresholds at every n-MAP MCMC realization. 
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To expedite this process but maintain the required level of 
precision, wc evenly down-sample long MCMC chains from 
the light curve fits until there are only 10® points remaining. 
Since the chains are evenly down-sampled, then the effective 
length of the chain is unaffected. 

For the n-MAP plots, plotting 10^ points on a figure 
is excessive both computationally and visually. Therefore, 
we down-sample any long n-MAP chains until we have l(f 
points remaining. Once again, the effective lengths are un- 
affected. 



4 HYPOTHETICAL SYNTHETIC SYSTEMS 

4.1 KOI-SOl: A Moderate-Eccentricity 
Triple-System 

In order to demonstrate and test the MAP techniques dis- 
cussed thus far, we will present two hypothetical analyses 
in this section. Additional control tests are also available in 
§A1, A2 & §A3. 

For our first example, wc consider a hypothetical three- 
planet system dubbed "KOI-SOl" for Kepler Object of In- 
terest Synthetic 01. We use the same identification as that 
used for Kepler Objects of Interest because we envision that 
real KOI targets will be the most obvious application of our 
technique in the near future. The three planets (KOI-SOl. 01, 
KOI-SOl. 02 & KOI-SOl. 03) are chosen to have orbital peri- 
ods of Pa = 13.9342 d, P2 = 25.13654 d and Pi = 44.86254 d 
around a Solar-star. These were selected to provide at least 
three transits for all three planets within the total time win- 
dow of the QO, Ql and Q2 Kepler data (127 d). We also 
deliberately avoid mean motion resonance to provide the 
plausible scenario that the planets follow a strictly linear 
ephemeris for the sake of simplicity (note that TTVs do not 
invalidate our technique and can be easily accounted for by 
allowing each transit to have an individual parameter for the 
time of transit minimum). The occultation depths are as- 
sumed to be negligible and the transit epoch for each planet 
is selected such that a) we obtain at least three transits for 
each planet b) no overlapping transits occur. 

The eccentricity parameters were selected such that the 
Hill stability criterion (Gladman 1993) was satisfied: 



+ 5.76 



2/3 



3/2 



(32) 



As a result of this criterion, the eccentricities are there- 
fore "moderate" and not large: 63 = 0.05, 62 — 0.15 and 
ei = 0.08. We also decided to enforce apsidal locking (Baty- 
gin et al. 2009) to try to provide potentially realistic sce- 
nario. Since the transit probability is highest for planets 
near to ~ 90° (Kane & von Braun 2008), we chose locking 
about this point: 0J3 = 92.23°, 0J2 = 89.42° and wi = 91.20°. 
Transit impact parameters were randomly generated to be 
63 = 0.276, 62 = 0.055 and 61 = 0.858 and planetary radii 
were arbitarily set to R3 = 0.45 P./, P2 = 1.05 P./ and 
Pi = 0.92 Pj. The properties of the KOI-SOl system are 
summarized in Table 2. 

Quadratic limb darkening coefficients for the star were 



generated assuming a Solar-like star and a Kurucz (2006)- 
style atmosphere, giving ui = 0.4277 and U2 = 0.2522. It 
should also be noted that the properties of all three planets 
satisfy the criteria for MAP to work, as described in §2.7 
i.e. the "If-cquation is accurate to better than 99% and the 
duration approximation is accurate to better than 99.9%. 

Synthetic data were generated to spam the 127 day win- 
dow of the QO, Ql and Q2 data of Kepler. We chose to 
use short-cadence (58.84876 s) data with Gaussian noise of 
250 ppm (consistent with typical Kepler noise for a V ~ 12 
star e.g. Kipping & Bakos 2011b) and random reference 
mean anomalies for the planets (although ensuring no mu- 
tual transits) yielding 186,393 synthetic photometry points. 

In a totally blind manner, one of us (DK) generated the 
synthetic data and blindly passed it onto the other three 
(WD, JJ & VM) who identified the number of planets, the 
orbital periods and fitted them using an MCMC routine 
coupled with the Mandel & Agol (2002) algorithm. We then 
followed the steps outlined in §3.1. 



4-1.1 Light curve fits 

The light curve fits for all three planets are shown in Fig- 
ure 4 and the derived posteriors are presented in Figure 5. 
Diagnostics on the mixing and convergence of these MCMC 
fits arc presented in Table 3, all of which indicate reliable 

results. 

As expected, the correct radii, transit epoch and 
orbital periods were found in the blind-search. The 
derived stellar densities, assuming a circular orbit, 

were found to be [pljf^'' = 1.47liE!;j;il g'^' cm^^^ 



circ\2/3 



1 -701 n+0 0052 2/3 „„-2 J ("^circx 

1.7019_o.oio9 g cm and (p,,3 ) 



2/3 



1.478to'o56 g^^^ cm ^' which clearly deviate significantly 
from both a common value and the actual stellar density 
of (p,)^^^ = 1.258 g^/^cm"^. Note that we here quote the 
median of each marginalized distribution as the best-value 
and the uncertainties come from the 34.15% quantiles either 
side of the median. This practice is continued for all results 
presented in this work. 



4.1.2 Results using a-MAP 

The analytic approximations from §2.6 may be used to pro- 
vide lower bounds on combinations of ei , 62 and 63 via Equa- 
tion 17, for which we find: (ei + 62) > 0.0779j:E;!;^^f, (true 
ei+e2 = 0.23), (63-^63) > 0.0752l'(;;;;j;^^ (true 62 + 63 = 0.20) 
and (63 + 61) > 0.0016t^:2?^^ (true 63 + ei = 0.13). As ex- 
pected, all of these values are consistent with the true num- 
bers. Further, the significance of each combination being > 
is given by 13.0-cr, 8.8-0" and O.l-a respectively, thus indicat- 
ing that the a-MAP method definitively shows that the sys- 
tem contains significant eccentricities. Given that (63 -I- 61) is 
consistent with zero, one may correctly assert that planet 2 
is most likely to contain the majority of the net eccentricity. 



4.1.3 Results using n-MAP 

Using the n-MAP algorithm described in §3.1, we explored 
the full 6-dimensional permitted parameter space with B = 
1.25 X 10® MCMC trials. Jump sizes were selected to be 1% 
for all terms (i.e. Aeu = 0.01, Acok = 0.01 x 27rrads). The 
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Table 2. Physical properties of the hypothetical KOI-SOl system. We assume Mp ^ M, during the calculation of the planetary 
semi-major axis, thereby removing the need to assign planetary masses. 



Object 


M 


R 


P 


b 


e 




KOI-S01.03 




5.1 R<s 


13.93 d 


0.28 


0.05 


92.2° 


KOI-S01.02 




11. 7i?© 


25.14 d 


0.05 


0.15 


89.4° 


KOI-SOl.Ol 




10.3 fl© 


44.86 d 


0.86 


0.08 


91.2° 


KOI-SOl 


1.00 Mq 


1.00 Rq 











Table 3. MCMC diagnostics of fits for KOI-SOl. Diagnostics presented here as discussed in §5.6. 



Planet 


# of Accepted 


Lowest efT. 


Parameter w/ 


Highest 


Parameter w/ 




MCMC trials 


length 


lowest eff. len. 


Geweke diag. 


highest Geweke diag. 


KOI-SOl.Ol 


1 X (1.25 X 10^) 


6247 


b 


0.028 


OOT_i 


KOI-S01.02 


1 X (2.50 X lO'^) 


1968 


b 


0.0051 


OOTo 


KOI-S01.03 


4 X (1.25 X lO'^) 


1258 


b 


0.032 


OOT+2 


n-MAP 


1 X (1.25 X 10*^) 


1156 


ei 


N/A 


N/A 




Figure 5. Marginalized posteriors from the light curve fits (stage 1) only, for KOI-SOl. In all figures, the green line marks the truth, 
with the numerical value provided in parentheses. Row 1: Posteriors for (p"^'^)'^^'' . Fits assume a uniform prior in 

(pf'kf^^- Notice 

how even though the eccentricities in the system are not large, the derived .stellar densities clearly indicate the presence of eccentricity 
in the system due to the non- overlapping posteriors. Row 2: Posteriors for the ratios of the (p"^°)^''^ terms. Row 3: Posteriors for the 
ratios of the minimum pair- combined eccentricities, computed using the approximations Eguation 11&13 (i.e. a-MAP). 
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Figure 4. Folded light curves of KOI- SO 1.01 (top), KOI- SO 1.02 
(middle) and KOI-S01.03 (bottom). The synthetic photometry is 
shown as gray circles, phase-binned with a bin size given by the 
number of transit epochs. The best-fit transit model is shown as 
a continuous solid line in each case. 



starting point for the chain was randomly generated until a 
point with < 10 was located. Mixing was checked for as 
described in §3.6 and the results are reported in Table 3. 

The results are shown in Figure 6, after down-sampling 
to one million trials. The joint posteriors clearly shows the 
minimum constraints on (ei + 62) and (e2 + 63) (the white 
regions in the corner), as derived using a-MAP. We note 
that the diamonds shown in the joint posteriors of Fig- 
ure 6 (which represent the true values), consistently lie 



in densely populated regions, supporting the validity of 
the MAP technique. The ID marginalized posteriors yield 
ei = 0.105t°:g?«, 62 = 0.138l°;g°4 and es = O.Uetngg. all 
consistent with the truth to within la. 

With the a-MAP results, it was shown that one may 
reasonably deduce that planet 2 is most likely to sustain 
the highest eccentricity in the system. n-MAP agrees with 
this conclusion since 62 has the highest median of all three 
marginalized posteriors. However, the significance of the or- 
bital eccentricities for all three planets appears marginal. 
The significances of e/s > for fe = 1, 2, 3 are l.A-a, 1.7-a and 
1.4-cr respectively, which a far cry from the 13.0-cr level de- 
tections which were achieved by a-MAP (on the same data). 
This suggests the following paradigm about the results from 
MAP (both a-MAP & n-MAP): MAP is more sensitive to 
pair- combined eccentricities than individual eccentricities. 

If we had not used n-MAP but assumed uniform pri- 
ors in efc, the probability that ei < 0.3 (which is a useful 
rough limit for a habitable world) would be 30% and the 
probability that ei > 0.3 would be 70%. Thus it would be 
~ 0.4 times more likely that the orbit was ei < 0.3 than 
otherwise. The same is of course true for 62 and 63. Using n- 
MAP these odds ratios become 2.4, 1.2 and 2.7 for ei < 0.3, 
62 < 0.3 and es < 0.3 respectively, demonstrating the extra 
information we have gained from using n-MAP. 

4.2 KOI-S02: 61-Virginis Analog System 

4.2.1 Setup 

As an additional test, we decided to look at a genuine mul- 
tiple planet system with well-characterized eccentricities. In 
order to satisfy this requirement and additionally locate a 
system with ^ 3 planets, we must draw upon planets found 
through radial velocity, rather than the transit technique. 
This is because RV planets have nmeh better orbital solu- 
tions than the few multiple systems found by Kepler so far. 
As the systems are RV planets and not transiting, we must 
generate synthetic photometry for them, in a similar way as 
to was done for KOI-SOl. 

From the 12 systems which satisfy our criteria at the 
time of writing (according to www.exoplanets.org), we se- 
lected the 61-Virginis system. 61-Virginis has orbital periods 
short enough to be detected within the first 18 months of 
operation of Kepler and posseses the highest orbital eccen- 
tricity components from such systems. Despite the eccentric- 
ities being the largest found from those available systems, 
all three planets in the 61-Virginis system satisfy the criteria 
for MAP to work, as described in §2.7 i.e. the ^-equation is 
accurate to better than 99% and the duration approximation 
is accurate to better than 99.9%. 

The relevant properties of the 61-Virginis are provided 
in Table 4 and taken from Vogt et al. (2010). We dub our 
hypothetical system as KOI-S02, to stress the fact that this 
is a hypothetical analysis and not a genuine study of 61- 
Virginis. 

In order to have guarantee that we have measured three 
transits of all three objects (assuming all three indeed tran- 
sit), we would require 4Pd days of continuous photome- 
try. To mimic this, we consider 1.5 years of Kepler short- 
cadence photometry. We estimated our noise based upon a 
simple calculation. The V = 11.4 star TrES-2 has been ob- 
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Table 4. Physical properties of the hypothetical KOI-S02 system, whose properties are assumed to be the same as that as the 61- Virginis 
system from Vogt et al. (2010). The radii and impact parameters have been arbitarily assigned, since the system is not known to have 
any transiting members. 



Object 


Object Analog 


M 


R 


P 




b 


e 




KOI-S02.03 


61-Vir b 


5.1 Af® 


i.6i?e 


4.215 


d 


0.15 


0.10 


110° 


KOI-S02.02 


61-Vir c 


10.5 M(B 


3.3 Ri^ 


38.02 


d 


0.40 


0.14 


340° 


KOI-S02.01 


61-Vir d 


22.9 Mm 


4.3 Re 


123.0 


d 


0.75 


0.35 


310° 



KOI-S02 61-Vir 0.94 M© 0.979 
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served in Kepler short-cadcncc mode to have an RMS noise 
of 237 ppm per minute (Kipping & Bakos 2011b). As 61- 
Virginis is V = 4.74, it would never have been observed 
by Kepler because it is too bright. Nevertheless, consider 
the star was at the bright hmit of Kepler^s range, namely 
F ~ 9, then the star would be at the floor-limit of RMS pre- 
cision. This should be around 25% lower than the RMS for 
a y = 11.4 star. Accordingly, we assigned an RMS precision 
of 178 ppm per minute for this synthetic data set. 

As the planets do not transit, wc had to assign planetary 
radii. We made the simple assumption that those planets of 
mass M > 10 (61-Vir c & d) were ice/gas giants sim- 
ilar in composition to Neptune. Wc accordingly computed 
their radius assuming the average bulk density was equal to 
PNeptunc ~ 1.628 gcm""^. For planets of mass (61-Vir h), we 
assumed the simple mass-radius scaling law of a terrestrial 
planet R ~ M^'^'' (Valencia et al. 2006) in Earth-mass and 
-radii units. 

Limb darkening was assumed to be the same as with 
KOI-SOl for simplicity. The transit impact parameters were 
arbitrarily chosen to be 6;, = 0.15, he = 0.40 and ha = 0.75. 



4-2.2 Light curve fits 

The light curve fits for all three planets are shown in Fig- 
ure 7 and the derived posteriors arc presented in Figure 8. 
Diagnostics on the mixing and convergence of these MCMC 
fits are presented in Table 5, all of which indicate reliable 
results. 

As expected, the correct radii, transit epoch and or- 
bital periods were recovered. The derived stellar densities, 
assuming a circular orbit, were found to be (p"'i^)^''^ = 
0.8021^03? g'^'cm-^ {pt^r^ = 1.035l°:S«tg'/'cm-2 
and {p''^}if'^ = l.SOgtfj-J^^g^/^'cnr^^ which clearly devi- 
ate significantly from both a common value and the actual 
stellar density of {p*f^^ = 1.260 g^/^ cm"^. 



4-2.3 Results using a-MAP 

The analytic approximations from §2.6 may be used to pro- 
vide lower bounds on combinations of ei , e2 and e-i via Equa- 
tion 17, for which wc find: (ei + cz) > 0.145to 054 (true 
ei -I- 62 = 0.49), (62-1-63) >0.215in^^ (true 62 -|- 63 = 0.24) 
and (63 + ei) > 0.425lS:g|| (true 63-1-61 = 0.45). As ex- 
pected, all of these values are consistent with the true num- 
bers. Further, the significance of each combination being > 
is given by 2.7-a, 2.4-a and 4.3-cr respectively, therefore the 
a-MAP method strongly indicates that the system contains 
non-zero eccentricities. 




-0.5 0.0 0.5 

Time, from transit m inim um [daysj 



1.0005 




^ 0.9990 



Ct^85 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 
Tinie ffiOJil tratisit minimiim [dayS] 




03997 



-0,3 -0.2 -0.1 0.0 0,1 0.2 0.3 
Thaisifsmn. transit minittHim [days] 



Figure 7. Folded light curves of KOI-S02.01 (top), KOI-S02.02 
(middle) and KOI-S02.03 (bottom). The synthetic photometry is 
shown as gray circles, phase-binned with a bin size given by the 
number of transit epochs. The best-fit transit model is shown as 
a continuous solid line in each case. 



4.2.4 Results using n-MAP . . . „ , . , . . / , 

mmimum constraints on all three pair-combmations (the 

Using the n-MAP algorithm described in §3.1, wc explored white regions in the corner), as derived using a-MAP. We 

the full 6-dimensional permitted parameter space with B = note that the diamonds shown in the joint posteriors of 

1.25 X 10® MCMC trials. Jump sizes were selected to be 1% Figure 9 (which represent the true values), consistently lie 

for all terms (i.e. Aej, = 0.01, Aw^ = 0.01 x 27rrads). The in densely populated regions, supporting the validity of 

starting point for the chain was randomly generated until a the MAP technique. The ID marginalized posteriors yield 

point with < 10 was located. Mixing was checked for as ei = Q.2Ata'S,^ ^2 = Q.2lt^^\\ and 63 = 0.321"; j?, all 

described in §3.6 and the results are reported in Table 5. consistent with the truth to within 1-cr except 63, which 

The results are shown in Figure 9, after down-sampling is marginally overestimated by 1.3-cr. Note that the mode 

to one million trials. The joint posteriors clearly shows the tends to overestimate the eccentricity because eccentricity is 
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Table 5. MCMC diagnostics of fits for K0I-S02. Diagnostics presented here as discussed in §5.6. 



Planet 


# of Accepted 


Lowest off. 


Parameter w/ 


Highest 


Parameter w/ 




MCMC trials 


length 


lowest eff. len. 


Geweke diag. 


highest Geweke diag. 


KOI-S02.01 


1 X (1.25 X 10^) 


4808 


b 


0.0048 


OOTo 


KOI-S02.02 


1 X (1.25 X 10^) 


2232 


b 


0.0074 


OOTo 


KOI-S02.03 


4 X (1.25 X 10^) 


1377 


b 


0.0027 


OOT+2 


n-MAP 


1 X (1.25 X 10^) 


1247 


£3 


N/A 


N/A 




Minimum {ci+t'2) Minimum (^i+i^-^) Minimum 3+^1) 

Figure 8. Marginalized posteriors from the light curve fits (stage 1) only, for KOI-S02. In all figures, the green line marks the truth, 
with the numerical value provided in parentheses. Row 1: Posteriors for (p'^^^l:)'^^^ ■ Fits assume a uniform prior in 
Posteriors for the ratios of the (p"^'^)^^^ terms. Row 3: Posteriors for the ratios of the minimum pair-combined eccentricities, computed 
using the approximatations Equation 11&13 (i.e. a-MAP). 



a positive definite quantity (Lucy & Sweeney 1971). These 
results, as with KOI-SOl (§4.1) support the validity of the 
MAP technique. 

In the earlier case of KOI-SOl (§4.1), Figure 6 showed 
that when the median of (e^ + ej)min posterior (derived us- 
ing a-MAP) was plotted along with the 2D-marginalized 
posteriors of ei versus (derived using n-MAP), there ex- 
isted excellent agreement on the minimum bound on the 
pair-combination. In contrast, the KOI-S02 simulation re- 
veals one interesting exception to this pattern, specifically 
for (ea -I- ei)min, as seen in Figure ??. Here, the median 
of the (ea + ei)niin posterior (derived using a-MAP) yields 
(ea -I- ei) 0.425 whereas visual inpsection of the n-MAP 
results suggests (ea -|-ei) 0.275. Including the a-MAP un- 



certainties reveals (ea -I- ei) 0.425to ggg and thus is 1.5-(t 
deviant from the n-MAP result. Whilst this is not statis- 
tically significant (and in fact both the a-MAP and the n- 
MAP limits are consistent with the truth of (ea-l-ei) = 0.45), 
this is still a large departure relative to the other cases, and 
may lead the reader to question the origin of this discrep- 
ancy. 

The reason for the discrepancy can be understood in 
terms of the approximations made in the original derivation 
of a-MAP. In §2.4, the derivation of a-MAP includes a step 
where we assume e^ <C 1, allowing us to execute a first- 
order series expansion. Therefore, one can see that the higher 
the eccentricity terms, the more a-MAP will depreciate as a 
reliable approximation. In the KOI-S02 case, the true values 
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of the eccentricity pair-combinations are (ei + 32) = 0.23, 
(62 + 63) = 0.20 and (63 + 61) = 0.45. Consequently, it is clear 
that (ea + ei) is exposed to the highest eccentricity and thus 
one should expect that a-MAP would be the least reliable for 
this case. Indeed, this is exactly what is seen, explaining the 
anomalous behaviour discussed in the previous paragraph. 
Despite the large pair combined eccentricity of 0.45, it is 
reassuring that a-MAP still performs quite well, landing to 
within 1.5-cr of the more sophisticated n-MAP result. 



5 DISCUSSION & CONCLUSIONS 
5.1 Recurring Patterns in MAP 

In the 2D posteriors shown in the example systems KOI- 
SOG, KOI-SOC, KOI-SOl and KOI-S02, it is possible to 
visually detect some recurring patterns. In this subsection, 
we will discuss the reason for these patterns. 

The most obvious pattern is that the plots tend to ex- 
hibit higher densities for ex — ey. This results in a diagonal 
high density region extending from the bottom-left to the 
top-right, getting narrower as it goes up. Why does this 
happen? 

Consider two planets on circular orbits. The light curve 
derived stellar densities will be identical meaning MAP will 
find that a system with two circular orbits to be a highly 
compatible solution. However, it is not the only compati- 
ble solution. Consider a MAP realization with ei = and 
e-z ~ 0.1. In this case, a rather broad range of UJ2 values can 
reproduce a stellar density which is still equivalent (or within 
error) to pS']". As 62 increases to higher eccentricities (i.e. as 
we move along one of the axes), the range of UJ2 values which 
can still reproduce such a result diminishes. Consequently, 
the integrated probability density over all 012 values is small 
and thus the region becomes a low-probability sector. 

Now consider MAP realizations where both ei and 62 
depart from zero. In this case, a situation where 62 = 0.9 and 
61 =0.1 is not grossly different from the case just described 
where 62 is very large and 61 is zero. Thus these regions 
tend to be low probability. However, if ei ~ 62 and both 
values are large then the range of wi and UJ2 values which 
can reproduce p'i^i ~ pl^f is broader. Thus the integrated 
probability density becomes larger too. This can be seen be 
simple inspection of Equation 8. For this reason, a common 
pattern seen in MAP is the tail extending along the = Cy 
axis. 



5.2 MAP in Combination with Other Observables 

Multibody Asterodensity Profiling (MAP) can provide con- 
straints on the orbital eccentricity of transiting planets in 
multiple planet systems without the need for radial ve- 
locities, occultations or transit timing variations. However, 
there may some cases where these observables are available. 
In such a case, MAP can be combined with these other ob- 
servables to further refine the constraints on the orbital ec- 
centricities. Our n-MAP technique utilizes an MCMC rou- 
tine to explore the parameter space of possible orbital eccen- 
tricities. This MCMC works by assigning a merit-function 
to each point and then proceeding via the Metropolis- 
Hastings rule. Using n-MAP alone, the merit function is: 



= Xmap (33) 

where Xmap is given in Equation 22. If additional ob- 
servables are available, such as TTVs or RVs then one may 
simply append the associated merit functions: 

= Xmap + Xrv + Xttv + -• (34) 

Once appended, the n-MAP routine is simply executed 
as before. In our next paper, we will provide an analysis of 
this technique on a real transiting planet system featuring 
both radial velocities and transit timing variations. 

5.3 Differences to the Moorhead et al. (2011) 
Technique 

A direct comparison to the method of Mil is not fair be- 
cause the techniques operate under different conditions and 
assumptions. Despite this, we will here outline a few differ- 
ences between the two methods. One advantage of the Mil 
method is that it works for all transiting planets whereas 
MAP requires > 1 transiting planet in a given system. One 
advantage of MAP over the Mil technique is that MAP is 
highly model independent'^, assuming only the planets orbit 
the same star whereas Mil rctjuirc a value for the stellar ra- 
dius determined via the more model sensitive route of stellar 
evolution. 

In principle, MAP exploits more information in the light 
curve than that of the Mil technique. Mil compare the ob- 
served transit duration to the maximum theoretical value 
for a circular orbit; in other words Mil make use of one 
metric for the duration. Typically, this metric is the transit 
duration defined as the time for the planet to move from its 
centre crossing the stellar limb to exitting under the same 
condition, T, since this is independent of the derived plane- 
tary radius. In contrast, MAP uses both T and the ingress 
duration, T12, to derive the light curve derived stellar den- 
sity assuming a circular orbit, pl"^''. Thus, it can be appreci- 
ated that MAP uses the same information as Mil plus some 
extra information. In the limit of ignoring this ingress infor- 
mation, the fundamental transit information used by both 
techniques would be identical. 

As a result of Mil negating the ingress duration infor- 
mation, the impact parameter is unresolved. For this reason, 
Mil adopt the conservative assumption that the limiting 
case is for b = 0. Ford ct al. (2008) alternatively discuss 
how a prior in h could be adopted by assuming an isotropic 
distribution in orbital inclination. In contrast, MAP takes 
the distribution in b from the data itself, essentially charac- 
terized by the ingress duration. Due to the 6 = conserva- 
tive assumption made by Mil. only transit durations longer 
than this limiting case can ever be detected. Unlike MAP, 
this limits the method to detecting eccentric planet transit- 
ing near apoccntre, the slowest part of the orbit, and also 
the least likely geometric configuration to detect a transiting 
planet in (Kane & von Braun 2008) (this is also pointed out 

^ There is a very weak model dependency in MAP via the 
adopted limb darkening law, but this can also be fitted for with 
sufficient signal to noise 
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Figure 9. Results of n-MAP fits for KOI-S02 using uniform priors in and uif;. Green diamonds/lines mark the truth. Blue lines 
mark the median of the lower limits found using analytic MAP (shown in Figure 8) Black denotes 0-1 a, red denotes l-2a, orange 
denotes 2-3a and yellow denotes > Sc. White regions were never vistted in any MCMC realizations and thus are highly improbable. e\ 
vs 63 plot has a displaced blue line due to the large uncertainty on this constraint, as evident in Figure 8. 



in Tingley & Sackett 2005). This bias, discussed in Mil, re- 
quires de-biasing any eccentricity statistics deduced and of 
course reduces the overall sample size since only a subset of 
eccentric planets are detected. For this reason, we anticipate 
MAP would provide a more powerful diagnostic of the statis- 
tics of eccentric planets, but of course is only applicable in 
multiple systems. 



5.4 Single-body Asterodensity Profiling (SAP) 

MAP makes no assumption about the properties of the par- 
ent star. For stars with poor characterization, this is an ad- 
vantage since the stellar properties are frequently subject to 
unknown systematic uncertainties. However, in some cases 
the stellar properties are well-characterized and this is in- 
formation which MAP ignores. An example of this is a star 
which has been studied with asteroseismology, leading to a 
highly precise determination of p, . 

Single-body Asterodensity Profiling (SAP) is the logi- 
cal extension of MAP which can include this information. 
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In this work, wc argue in favour of not using SAP due to 
the frequently unreliable stellar parameters and the strong 
model dependency of those derived results, which is in stark 
contrast to the MAP technique. However, for sake of com- 
pletion we will discuss here a possible implementation of 
SAP. 

In each MCMC trial of the n-MAP algorithm, we create 
a {ek,i^k} trial vector. This may be used to construct a 
{^Pfc} vector. For a single planet, this would be simply be a 
1-dimonsional vector: {^i}. This value may be used to infer 
the true stellar density based upon the derived value of p^'f 
via: 



„SAP „circ /,x, 
P*,l = P*,l / Vl 



(35) 



For each MCMC trial, we can use the same value of Pf"i, 
namely the median of the light curve derived posterior. We 
now simply define xIap '^s the number of standard devia- 
tions between this trial value of p^^i^ and the empirically 
determined value of p* , say p! 



seismology . 



XSAP.fe — 



lP*,k — P* 



seismology ] 2 



'yiplf)? + [^{pT 



seismology 



(36) 



where we have replaced the 1 subscript (for planet "1" ) 
with k to make it a general equation. Clearly for n planets 
we have n contributions to the total function from SAP: 



2 _ \ " 2 

XSAP — XSAP,fc 



(37) 



a synthetic system with one planet with e = 0.15 (KOT 
SOl, see §4.1), Q0-Q2 Kepler data can infer a significant 
eccentricity at the 13-cr level. 

To determine the individual eccentricities require the 
use of numerical methods, or n-MAP. n-MAP is shown 
to recover the same constraints for the pair-combined ec- 
centricities as a^MAP does, but additionally provides indi- 
vidual constraints as well as a fuller picture of the inter- 
relationships between the various eccentricity terms. How- 
ever, wc find MAP is more sensitive to the minimum pair- 
combined eccentricities than individual terms. Nevertheless, 
an empricial determination of the posterior for each term is 
always derivable. The highly model independent nature of 
MAP means that these posteriors will narrow ad-infinitum 
as more data and signal-to-noise is accumulated. Further, 
the application of MAP to dozens of systems raises the po- 
tential of discerning truly unbiased statistics on the eccen- 
tricity distribution of planets in multiple systems. 

Ultimately, MAP has the potential to characterize the 
eccentricity of the first truly habitable Earth-like planet, 
which could be found by Kepler. In such a case, the Earth- 
mass planet will likely be too challenging to detect with ra- 
dial velocities but MAP will be an ever-present tool provided 
the photometric time series is of reasonable quality and the 
system has more than one transiting planet (true of ~50% 
of all transiting planet candidates found by Kepler) . In con- 
clusion, MAP offers a method to diagnose both interesting 
dynamical systems and potentially interesting astrobiologi- 
cal targets too. 

The n-MAP algorithm is available as a Fortran 90 script 
upon request. 



This may be combined with the MAP, RV or TTV merit 
functions as desired to produce finer constraints on the or- 
bital eccentricity. 



5.5 Overview 

In this work, we have presented a new method to photo- 
metrically constrain the orbital eccentricities of transiting 
planets. The method is only applicable to multiple transit- 
ing planet systems and relies on the key assumption that 
all of the transiting planets orbit the same star. The new 
method works by comparing the light curve derived stellar 
density between each planet and thus is dubbed "Multibody 
Asterodensity Profiling" (MAP). MAP requires no prior in- 
formation on the star's properties and thus is highly robust 
against systematic uncertainties. 

MAP constitutes a new observable for which the likeli- 
hood of a given orbital configuration can be computed in a 
X^-sense. Thus, MAP be be combined with other pieces of 
information about the eccentricity of the system e.g. tran- 
sit timing variations, radial velocities, occultations. In this 
work, we have adopted uniform priors for the orbital eccen- 
tricities but more realistic priors based upon dynamics or 
planet formation can also be invoked (see §3.3). 

In its simplest form, MAP can be applied by employing 
some simple analytic expressions (a-MAP) to deduce the 
combined eccentricities of two planets (e.g. Equation 11). In 
this sense, MAP exhibits impressive sensitivity to a system 
with even a moderate to low eccentricity. For example, for 
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APPENDIX A: FURTHER HYPOTHETICAL 
SYNTHETIC SYSTEMS 

In this appendix, wo will present several additional test sys- 
tems on which we implemented the MAP technique. The 
systems are control cases uses to ensure we retrieve the cor- 
rect null results. 



Al KOI-SOP: Triple System Fit Driven by 
Uniform Priors Only 

A 1.1 Setup 

We repeated our analysis of KOI-SOl (see §4.1 for details) 
but in the second MCMC chain instructed the algorithm to 
accept all jumps, regardless of the likelihood. This instruc- 
tion causes the algorithm to produce a result which purely 
reproduces the initial priors and thus is a useful of way of 
visualizing what an n-MAP result appears like in the ab- 
sence of any observational constraints. We dub this system 
KOI-SOP (for "priors"). The resulting n-MAP posteriors are 
shown in Figure Al, which clearly reproduce the expected 
behaviour of uniform priors. 

A2 KOI:S0G: Triple System with Normally 
Distributed {pf'^f^'^ 

A2.1 Setup 

For our second control system, we generated three inde- 
pendent distributions for {p"^kY^^ of 10* points each. Each 
(pS'fe^)'^^^ realization is generated with random numbers se- 
lected from a normal distribution centred about n = 1 with 
standard deviation a = 0.01. This allows us to see the iso- 
lated behaviour of MAP without any concern for potential 
biases from the light curve fitting stage. The resulting pos- 
teriors arc shown in Figure A2. 

The n-MAP MCMC stage stopped when 3.75 x 10*^ tri- 
als had been accepted, using 1% jump sizes for all six free 
parameters (ei, wi, eg, W2, es, W3). Mixing was checked by 
evaluating the effective lengths of all free parameters from 
the non-burnt trials, and we found the lowest effective length 
was 1098 for parameter ei. In order to keep only 1.0 x 10^ 
points in the figures, we evenly down-sampled the chain to 
produce the figures. Table Al summarizes the MCMC diag- 
nostics. 

Since all values were selected to be equal, we expect the 
system to be compatible with a triple-circular orbit system. 
The n-MAP results indeed agree with this conclusion, as 
shown in Figure A3. The bulk of the black points (< l-tj) 
land close to circular orbit solutions and the marginalized 
posteriors of eu reflect the strong preference towards a low- 
eccentricity system. 

Both the a-MAP and n-MAP results support the con- 
clusion of a near-circular triple system. The 90% upper lim- 
its on the eccentricity were found to be ei < 0.31, 62 < 0.28 
and 63 < 0.33. 

If we had not used n-MAP but assumed uniform priors 
in e^,, the probability that ei < 0.3 (which is a useful rough 
limit for a habitable world) would be 30% and the proba- 
bility that ei > 0.3 would be 70%. Therefore it would be 
~ 0.4 times more likely that the orbit was ei < 0.3 than 
otherwise. The same is of course true for 62 and 63. Using n- 
MAP these odds ratios become 8.4, 10.4 and 7.8 for ei < 0.3, 
62 < 0.3 and 63 < 0.3 respectively, demonstrating the extra 
information we have gained from using n-MAP. 

A3 KOI-SOC: A Perfectly Circular Triple-System 

As a final test, we considered a system of three transiting 
planets, all on circular orbits, denoted KOI-SOC. This sys- 
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Table Al. MCMC diagnostics of fits for KOI-SOG. 



Planet 


# of Accepted 
MCMC trials 


Lowest eff. 
length 


Parameter w/ 
lowest eff. len. 


Highest 
Geweke diag. 


Parameter w/ 
highest Geweke diag. 


KOI-SOG.Ol 


1 X (1.00 X 10^) 


N/A 


N/A 


N/A 


N/A 


KOI-S0G.02 


1 X (1.00 X 10^) 


N/A 


N/A 


N/A 


N/A 


KOI-S0G.03 


1 X (1.00 X 10^) 


N/A 


N/A 


N/A 


N/A 


n-MAP 


1 X (3.25 X 10^*) 


2978 


ei 


N/A 


N/A 



22 Kipping et al. 




Minimum ((^i+t'i) Mmimum (ei+^i) Minimum (e^+ei) 



Figure A2. Marginalized posteriors based upon synthetic (p"'^P)^''^ distributions, chosen to be perfectly Gaussian (KOI-SOG). In all 
figures, the green line vertical line marks the truth, with the numerical value provided in parentheses. Row 1: Posteriors for (p"';f )^'^^- 

These have been synthetically generated to be perfectly Gaussian. Row 2: Posteriors for the ratios of the (p"^'^)^''^ terms. Row 3: 
Posteriors for the ratios of the minimum pair- combined eccentricities, computed using the approximate expressions Equation 11&13. 



tem is identical to KOI-SOl (see §4.1) except tiiat efc = for 
all k. The system is generated, noised, re-fitted and treated 
with a-MAP and n-MAP with precisely the same methodol- 
ogy used on KOl-SOl. The MCMC diagnostics are presented 
in Table A2, which indicate excellent mixing and conver- 
gence in all cases. 

As expected, the correct radii, transit epoch and 
orbital periods were easily found in the blind-search. 
The derived stellar densities, assuming a circular or- 
bit, were found to be (pl'ff^^ = 1.2551° °^^ g^/^ cm'^, 
ipiyf^' = 1.2583lHgig'/'cm-2 and ' (p?;^)^/^ = 
1.252l!^:2:^| g^/^ cm-2 (truth is 1.411 gcm"^). 

The n-MAP MCMC stage stopped when B = 2.5 x lO'^ 
trials had been accepted, using 1% jump sizes for all six 
free parameters (ei, uji, 62, tJ2, 63, uus). Mixing was checked 
by evaluating the effective lengths of all free parameters 
from the non-burnt trials, and we found the lowest effec- 
tive length was 1116 for parameter ei. In order to keep only 
1.0 X 10'' points in the figures, the chain was again evenly 
down-sampled. 

The results, shown in Figure A4&A5, appear similar to 
those of KOI-SOG (§A2), which is to be expected given the 
fact both cases are consistent with a triply-circular system. 
The results again reflect a range of solutions consistent with 
circular orbits. 



As expected, the MAP results suggest that all three- 
planets are consistent with a circular orbit. The 90% upper 
limits on the eccentricity are ei < 0.38, 62 < 0.35 and 63 < 
0.41. 

If we had not used n-MAP but assumed uniform pri- 
ors in efc, the probability that ei < 0.3 (which is a useful 
rough limit for a habitable world) would be 30% and the 
probability that ei > 0.3 would be 70%. Therfore it would 
be ~ 0.4 times more likely that the orbit was ei < 0.3 than 
otherwise. The same is of course true for 62 and 63. Using n- 
MAP these odds ratios become 5.8, 6.8 and 5.0 for ei < 0.3, 
62 < 0.3 and 63 < 0.3 respectively, demonstrating the extra 
information we have gained from using n-MAP. 

This paper has been typeset from a T^jX/ file prepared 

by the author. 
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(0.00: 1.3 tr) 
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Figure A3. Results of the numerical MAP fits for KOI-SOG, using uniform priors in ej^ and uj^. In all figures, the green lines and 
diamonds mark the truth. Black denotes 0-1 a, red denotes l-2a, orange denotes 2,-3a and yellow denotes > 3cr. Blue lines mark the 
median of the lower limits found using analytic MAP (shown in Figure A2), White regions were never vistted in any MCMC realizations 
and thus are highly improbable. 



Table A2. MCMC diagnostics of fits for KOI-SOC. 



Planet 


# of Accepted 


Lowest eff. 


Parameter w/ 


Highest 


Parameter w/ 




MCMC trials 


length 


lowest eff. len. 


Geweke diag. 


highest Geweke diag. 


KOI-SOC. 01 


1 X (1.25 X W^) 


7812 


b 


0.0034 


T 


KOI-SOC. 02 


1 X (1.25 X 10^) 


12499 


b 


0.0017 


T 


KOI-SOC. 03 


10 X (1.25 X 10^) 


1083 


b 


0.016 


OOT_3 


n-MAP 


1 X (2.50 X lO**) 


2442 


62 


N/A 


N/A 
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Figure A4. Marginalized posteriors from light curve fits (stage 1) only, for KOI-SOC. In all figures, the green line vertical line marks 
the truth, with the numerical value provided in brackets. Row 1: Posteriors for fitted (p"^"^)^/^. Row 2: Posteriors for the ratios of 
the (p"^'^)^''^ terms. Row 3: Posteriors for the ratios of the minimum pair-combined eccentricities, computed using the approximate 
expressions Equation 11&13. 
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ei = O.076_o,o59 
(0.00: 1.3 a-) 
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Figure A5. Results of the numerical MAP fits for KOI-SOC, using uniform priors in e^ and u)f^. In all figures, the green lines and 
diamonds mark the truth. Blue lines mark the median of the lower limits found using analytic MAP (shown in Figure A4)- Black denotes 
0-1 a, red denotes 1-2(7, orange denotes 2-3(7 and yellow denotes > 3 a. White regions were never vistted in any MCMC realizations 
and thus are highly improbable. 



